Vortex energy and vortex bending for a rotating Bose-Einstein condensate 
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For a Bose-Einstein condensate placed in a rotating trap, we give a simplified expression of 
the Gross-Pitaevskii energy in the Thomas Fermi regime, which only depends on the number and 
shape of the vortex lines. Then we check numerically that when there is one vortex line, our 
simplified expression leads to solutions with a bent vortex for a range of rotationnal velocities and 
trap parameters which are consistent with the experiments. 
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INTRODUCTION 



Since the experimental achievement of Bose-Einstein 
condensates in confined alkali-metal gases in 1995, there 
has been a huge experimental and theoretical interest in 
these systems @, | | | | | |, | | 0. The study 
of vortices is one of the key issues. Two different groups 
have obtained vortices experimentally, the JILA group 
J§ and the ENS group (j) [§. In the ENS experiment, 
a laser beam is imposed on the magnetic trap holding 
the atoms to create a harmonic anisotropic rotating po- 
tential. For sufficiently large angular velocities, vortices 
are detected in the system. Experimentally, the ENS 
group |8| has observed that when the vortex is nucleated, 
the contrast is not 100% which means that the vortex 
line is not straight but bending. Numerical computa- 
tions solving the Gross-Pitaevskii equation |n], Q have 
shown that there is a range of velocities for which the 
vortex line is indeed bending. The aim of this paper is 
to justify these observations theoretically in the Thomas 
Fermi regime. We define an asymptotic parameter which 
is small in the Thomas Fermi regime and approximate 
the Gross-Pitaevskii energy to obtain a simpler form of 
the energy which only depends on the shape of the vortex 
lines. Then we check numerically that our characteriza- 
tion leads to solution with a bent vortex for a range of 
values of the rotationnal velocity which are consistent 
with the ones obtained numerically [ [Tl|| . Let us point 
out that Svidzinsky and Fetter |f| have studied the dy- 
namics of a vortex line depending on its curvature. For a 
vortex velocity equal to 0, the equation obtained in |l3| is 
the same as the equation corresponding to the minimum 
of our approximate energy, though the formulation in |D| 
was not derived from energy considerations. Moreover, 
their analysis is only valid for a single vortex line. 

The Gross-Pitaevskii energy provides a very good de- 
scription of Bose-Einstein condensates: it is assumed that 
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the N particles of the gas are condensed in the same 
state for which the wave function <f> minimizes the Gross- 
Pitaevskii energy. In the ENS experiment, a laser is ap- 
plied to the trap which makes it rotate. By introducing 
a rotating frame at the angular velocity = fle z , the 
trapping potential becomes time independent, and the 
wave function <fr minimizes the energy 
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under the constraint J \cj)\ 2 — 1. Here, for any complex 
quantities u and v and their complex conjugates u and 
V, (u,v) = (uv + uv)/2. 

We want to nondimensionnalize the energy in order 
to get a parameter which is small in the Thomas-Fermi 
regime. This framework of study has been developed by 
one of the authors in juj, except that [H was a two 
dimensionnal study for a condensate confined in the z 
axis. We define the characteristic length d — (h/muJx) 1 ^ 2 



and assume u,. 



acu x , cu z 



(3lu x . We set 
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where g^o = iTrh 2 a/m. For numerical applications, we 
going to use the experimental values of the ENS group 
m = 1.445. 10- 26 kg, a = 5.8.10~ n m, N = 1.4.10 5 
and uj x = 1094s" 1 with a = 1.06, = 0.067. We find 
that e = 0.0174, thus, e is small, which will be our asymp- 
totic regime. We re-scale the distance by R — d/^/e 
and define u(r) = i? 3 / 2 </>(x) where x = Rr and we set 
fl = Cl/eu x . The velocity fl is chosen such that O < 1/e, 
that is the trapping potential is stronger than the inertial 
potential. The energy can be rewritten as: 



E 3D (u) = J ^|Vu| 2 + ft ■ (iu, Vm x r) 



(2) 
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Due to the constraint J \u\ 2 = 1, we can add to E^d any 
multiple of J \u\ 2 so that it is equivalent to minimize 
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fl ■ (iu,Vu x r) + -r^\u\ - —zp T r(r)\u 
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where p TP (r) = po — (^ 2 + ot 2 y 2 + f] 2 z 2 ) for some con- 
stant p to be determined. Let V be the ellipse {p TF > 
0} = {x 2 + a 2 y 2 + (3 2 z 2 < p }. We impose the following 
constraint on p TF : 



p TF (r) = 1. 



(3) 



Indeed, as e tends to 0, the minimizer will satisfy that \u\ 2 
will be close to p TF so that the constraint will be satisfied 
automatically by u if we impose (^). In other words, p TF 
is the Thomas Fermi approximation of u. Equation (j^) 
leads to 
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>o /2 = 15a/3/87r. 
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To study the problem analytically, it is reasonable to min- 
imize the energy over the domain T> with zero boundary 
data for u. Indeed, when p TF < 0, the energy is convex so 
that the minimizer u goes to zero exponentially at infin- 
ity (see the numerical observation in [|| and the analysis 
on the behavior near the boundary of T> as well as the 
decay at infinity of the order parameter in ) . We 

consider the problem 



mm E e {u) subject to u G H^(V), / \u\ 2 = 1 

Jv 
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where 



E £ (u) = ( ^|Vw| 2 + ft ■ (iu, Vu x r) 
Jv 2 

+ i ^(p TF (r) - |u| 2 ) 2 . 

Note that a critical point u of E e is a solution of 

1 



(5) 



Au + 2z(fi x r).V« : 



MPi 



p £ u in 2?, (6) 



with w = on &D and p £ is the Lagrange multiplier. The 
specific choice of po in (^|) will imply that the term p e u 
is negligible in front of p TF u/e 2 . 

We have set the framework of study of our energy. In 
Section 2, we will make an asymptotic development of the 
energy taking into account that e is small (but |log e| 
is not big). Then in Section 3, we will check that our 
approximate energy yields a solution which is consistent 
with the numerical and experimental observations. 



II. ASYMPTOTIC DEVELOPMENT OF THE 
ENERGY 

Our aim is to decouple the energy into 3 terms: a 
part coming from the solution without vortices, a vortex 
contribution and a term due to rotation. 



A. The solution without vortices 

Firstly, we are interested in solutions without vortices, 
that is u has no zero in the interior of V. Thus we con- 
sider functions of the form -q = fe lS , where r\ is in Hq(T>) 
and / is real and has no zero in the interior of V. We 
consider first minimizing E s over such functions without 
imposing the constraint that the L 2 norm is 1, that is, / 
and S minimize 

£ E (f,S) = J^l\Vf\ 2 + -L(p TF -f 2 ) 2 

+ \ J / 2 |V5-^xr| 2 -/W, (7) 
where r = xe x + ye y . We have f e = on dT> and 
-A/ £ +/V# £ (VS £ -2nxr) - \f e (p^-f!) in V, (8) 



div(/ 2 (V5 £ -nxr)) = 0. 



(9) 



Equation ([)]) implies that there exists £ e in H 2 {T>) fl 
#o(X>) such that 



/ 2 (V5 £ -f2xr) = f7curl £ E 



(10) 



So £ e is the unique solution of 



curl (-J curl f e ) 

/ e 



-2 in V, Ce = on &D. (11) 



In the special case where the cross section of V is a disc, 
the minimum of (|^) is reached for ViS = but this is 
not the case if the cross section is an ellipse and there is 
a non trivial solution of (g). As e tends to 0, since the 
ellipticity of the cross-section is small, f 2 tends to p TF in 
every compact subset of T> and the function £ e given by 
© or |Ll]) tends to the unique solution £ of 

curl (—curl f) = -2 in V, £ = on dV. (12) 



One can easily get that 



-p 2 TF {x,y)/{2 + 2a 2 )e 2 



(13) 



Using (|l0|), we can define So, the limit of S E , to be the 
solution of p TF (ySo — f2 x r) = fl curl £ with zero value 
at the origin. We have So = Cflxy with C — (a 2 — 
I) /(a 2 + 1). We see that Sq vanishes when a = 1 that 
is when the cross-section is a disc. This computation is 
consistent with the one in H, though it is derived in a 
different way. 

B. Decoupling the energy 

Let r] e = j e e lS ' be the vortex free minimizer of E e 
discussed previously without imposing the constraint on 
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the norm of u. Let u e be a configuration that will mini- 
mize E e and let v £ — u e lrj e . Since r\ e satisfies the Gross 
Pitaevskii equation we have 

/ p (M 2 - - - fl) + |V/ e e«* | 2 

-2/|(V5 £ n x r)) =0. 

Using this identity, one can get that the energy E e (u e ) 
decouples as follows 



E s (u s ) = E e (r) s ) + G Ve (v £ ) + (14) 



where 



and 



The first term in the energy is independent of the solution 
it e , so we have to compute the next two and find for 
which configuration u e the minimum is achieved. We 
assume that the solution u e has a vortex line along 7 
and we call S 7 the dirac measure along the line. Our 
aim is to estimate the energy of u £ depending on 7. A 
first approximation consists in writing that |77 e | 2 tends to 
p TF when e is small so that we can approximate G Ve by 
Gvp TF = G e and I Ve by 1^ = I e . 



C. Estimate of G e (v s ) 



We want to estimate 



G e {v e ) = 



D 



2 \VvA 2 
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The mathematical techniques to approximate G e have 
been introduced by 1 17 in the 2 dimcnsionnal case and 
by Jl8[ in dimension 3, when e is very small. The prob- 
lem here is that e — 0.0174 so that |log e| is not large and 
there will be additional terms in the asymptotic expan- 
sion. For a minimizing configuration, one can prove that 
v e is tending to 1 everywhere except close to the vortex 
line 7. We define 



T Xe = {xeV s.t. dist(x,7) < Ae}, 



(15) 



and assume that Ae is small, A being a nondimensionnal 
parameter to be fixed later on. Then we split G e into 
two integrals: one in T\ £ and the other in T> \ T\ e . 



1. Estimate near the vortex core 

We are going to estimate G E in T Ae . At each point -f(t) 
of 7, we define II -1 (7^)) to be the plane orthogonal to 7 



at "f(t). Since Ae is small, we assume that p TF is constant 
in n _1 (7(i)) n T\ e and we call the value p t — p TF (^(t)). 
We want to compute 



/ ^ TF |v, e | 2 + gf(i-M 2 ) 2 



Pt 



, , |v, £ | 2 + |,(i-M 2 ) 2 . 

1 ./n-i( 7 (f))nT ie Z£ 



Since u £ is a minimizing configuration of E e , after scaling 
by r^fpl/e, we find that v £ is very close to u\{r ^fpl / 1 e) , 
where ui(r, 9) = fi(r)e l6 is the solution with a single zero 
at the origin of 



Au + u(l - M 2 ) = in IR 2 . 



Thus, 



n-!( 7 (t))nT Ae 



|v« £ | 2 + ^(i-kl 2 ) 2 



= I |v Ul | 2 + i(i-M 2 ) 2 

~ + 27rlog(A % /pT), 
where 




• = i Si + 5(1 - ft? + 

H 2 1 



ft -I 



fl 



The last approximation is good if X^fpl is large (in fact 
bigger than 3 is enough). The existence of A is justified 
by the fact that y/p TF /e is much bigger than 1, except 
very close to the boundary. 

The final estimate of this section is 

P TF (y + irlog(Xyfp~^j)dl (16) 



2. Estimate away from the vortex core 

We are going to estimate G e in 2? \ T\ e . In this region 
v e ~ 1 , so that only the kinetic energy of the phase has 
a contribution. 



/ ip TF |v« £ | 2 + ff(i-M 2 ) 

- / ip TF |V0 £ | 2 , 

JV\T X , 1 



IV\T X 

where 4> e is the phase of v e . Of course, 4> e is not defined 
everywhere. But let ip be such that div ip = and 

curl ip = p TF V(j). 
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Then ip is the unique solution of 



1 



curl curl ip = 2tt<L, V = on &D, (17) 

V p TF / 

where <5 7 is the vectorial dirac measure along 7 and 
/ iU P |V<M 2 = / ^-|curl^| 2 

ip ■ V(f) e x v 



IT>\T; 
1 

2 



where v is the outward unit normal. We will see that ip 
is almost constant at a distance Ae from 7 and we call 
this value ip\ e {^). Since the vortex line has a winding 
number 2ir, 



1 



V\T X 



-p TV \Vct> E 
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We have to compute ip on dT\ s . The computation is 
inspired by the paper of Svidzinsky and Fetter |l|. It 
follows from ( |l7| ) that ^ satisfies 



-A0- 



Ptf 



X CUrl ?/> = 27T/9 TF 5 7 . 



Let Xq <E 7. We denote by e3 = 7(^0) an d (ei, 62,63) 
an orthogonal base in local coordinates. Then ip has 
coordinates ipi in e,; and the variations of -03 are the only 
ones of influence in the equation for ip, since we want to 
compute ip ■ dl. We approximate the equation for ip by 



AV-3 



V03 



27rp TF o 7 , 



(18) 



where p TF (x 1 , rr 2 ) = j o TF (x 1 , x 2 , ccq). Let 3 = ips/^fp\ 
Then it follows from (Jig) that S satisfies 



where 



AS + /xH = 2ir x /p TF 5 1 



^ A 1 , . 1 

/i = v /p TF A-== = ^fp^l\±—= 

VPtf v>i 



(19) 



(20) 



Here Aj_ is the laplacian in the plane perpendicular to 
63 = 7(^0). If the cross-section is a disc one can compute 
p. We denote by 9 the angle of e 3 that is e 3 = cos 9e r + 
sin(9e 2 and (r, z) are the coordinates of xq in the original 
frame. Then 



(1 



1 cos 2 6 3(r sin 6 — f3 2 z cos ( 



Ptf 



Ptf 



Note that p > 0. In fact our numerical computations 
even yield p > 7. Our aim is now to give an approximate 
expression for S. We locally approximate the curve 7 
near the point xq by the parabola x = kz 2 /2, where k is 
the curvature of 7 at x$. This is where we use the same 
ideas as in p3|. Note that in our approximations, we are 



only taking into account the shape of 7 close to xq. The 
justification for this relies on the fact that p > 7 as our 
numerics show. Indeed if we solve 

~AX + pX = f 

where / is supported at a distance d of xq. Then using 
the Green function, we find that 



\X\ < 



4:irp 3 d 

In particular, for d = 0.1, this gives an error less than 
10~ 3 . This is to be compared to the Euler constant and 
our approximation is reasonnable. We rewrite ( fl9| ) in 
local coordinates to get 



-As + kd xl E + pE = 277y/p TF ([x )5 e3 , 

where S e3 is the dirac mass supported along the line 63. 
Thus 



-A(e-^-»)+((-) (e~^"s) = 2Tr v / p TF (x )S e3 . 

The solution of this equation is 



Vft^N^o ( y P + — dist(x, 7) 

where Kq is a modified Bessel function. In particular, 
K Q {x) ~ -log(e c °x/2) for small x where C ~ 0.577 is 
the Euler constant. Hence, we deduce 



/ eP° I k? \ 
ip(x) ~ -p TF log(^ — W p + — dist(x, 7) J 7. (23) 

Thus we conclude by the estimate for G £ (v £ ) u\V\T\ £ 



G E {Ve)\V\T Xz ^ -7T / Pxrlog 

7 



(tV^ + t A£ ) dL (24) 



D. Estimate of I s (v e ) 

We want to estimate 

I e (v e )= I p TF (\/S E -flxr)-(iv s ,Vv E ). (25) 

Recall that the unique solution of ( fTT| ) satisfies p TF (V5, 
O x r) = curl £ e . Hence we integrate by part in (|2 
to get 

le(w e ) = O / £ e • Curl (w e , Vu e ). 

Jt> 

Let </> £ be the phase of v E . Since v E is tending to one 
everywhere except on the vortex line, then (iv E ,\7v E ) ~ 
V</> £ , hence we can approximate curl [iv E , \7v E ) by 2n5j. 
We use the value of £ given by ( |l3| ) and the fact that 
j(t) ■ e z = dz to get 



Ie(Ve) 



(1+a 2 ) 



Ar F dZ- 



(26) 
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E. Final estimate for the energy 



We use (y|)-(|l6|)-(gj)-(g6|) to derive the energy of a 
solution with a vortex line. Indeed the energy of any 
solution minus the energy of a solution without vortex is 
roughly the vortex contribution in the sense: 



E E (u e ) - E e {r) e ) ~ £, 



We find that the vortex contribution £L is 



(27) 



Ptf ( y + 



7Tl0g( 



Ptf 



f// 



(1 + a 2 ) 



2 j 
2 \ / Ptf 



(28) 



Hence if the right-hand-side of (|2^) is negative, it means 
that it is energetically favorable to have vortices. Note 
that in the first integral of £ 7 , we have dl = |7(z)|dz 
whereas in the second one, we have dz. 

If the vortex line is straight, our computation yields 



3/2 

Po 



V2 



5f 



Co) 



2tt. 



/-10 

-7T 

V 9 



-log2 



n- 



3 

87rp 



logpo 



15(1 + a 2 ) 



(29) 



This gives a critical angular velocity f^i for which a 
straight vortex has a lower energy than a vortex free so- 
lution. With our experimental data, it yields f^i ~ 22.45, 
that is Cl\/u> x ~ 0.39. We are going to see in the numeri- 
cal section that for Q < Cli, a bent vortex has a negative 
energy. 



F. Case of several vortices 

Let us assume that the solution u s has n vortices along 
the lines 74, 1 < i < n. We want to estimate the energy 
in this case. For each 7$, we define Tj a £ as i n (|15| ). 

One can check that the estimates ( p6| ) and (|16|), respec- 
tively for I € (v 6 ) and for G E (u e ) close to each vortex core, 
are unchanged if the integral along 7 is replaced by the 
sum of the integrals along 7,. The only difference is for 
the estimate away from the vortex cores where we have 
to take into account the interaction between the vortex 
lines. Let us denote T> n =T>\ \JjTj t \ e . We still have 



G e (v e ) 



1 



x> 2p TF 



■|curl ip\ 2 , 



(30) 



where ip = y^. ipj and ipi solves ( p"7j ) with 7; instead of 7. 
Thus, we need to estimate 



E 



V, 



2 P i 



I curl ^| 2 +E 



2p TF 



I curl V'fcl • |curl ^>j|. 

(31) 



The first integral is estimated as in section C.2 by 

0C0 / £*2 



\ f /e / k 2 \ 

Y,- % j P T Flog(— ^Jh + —XeJ dl. (32) 

As for the second integral in ([jj]), we integrate it by part 
to get 



(33) 



The computation of ipk(%) from section C.2 is still valid 



and we have ipk{x) — — PtfKo( y A* + ^-dist(x, Jk))- This 
yields the contribution of n vortex lines (to be compared 
with (E8[) for 1 vortex) 



£n= E/ ^(y+Trlog 

Ptf rf 2 



'74 
^7T 



Ptf 



Ee Co W M + f 



(1 + a 2 ) 



(34) 



y p TF K (^j fx + ^dist(x,jk))dl, (35) 



where _Ko is a modified Bessel function. Note that the 
curves are going to interact only in the region where they 
are close to one another. 



III. COMPARISON WITH THE NUMERICS 

We are interested in the shape of the vortex line that 
minimizes ( p8[ ) according to the value of fl. We write 
(7(t) = r{t),z{t)) and we assume that the vortex line is 
in the plane (y, z). We will denote p TF (t) — p — a 2 r 2 (t) — 
(3 2 z 2 (t) and we define 



lp -a 2 r 2 {t)~ f3 2 z 2 {t) 

£.2 



Since (|28j) does not depend on the parametrization 7(t), 
we choose a special parametrization on the curve such 
that 



C 2 (t)p 2 TF (t)^(t) + z 2 (t)) = l 
Then it is equivalent to minimize 

. .0 , . . .> , , , . Qir 



C 2 (t)p 2 TF {t)(r 2 {t) + z 2 (t))dt- 



(1 + a 2 ) 



(36) 



p 2 TF (t)z(t)dt 



(37) 

under the constraint (|36|). In our computations below, 
we will proceed to a minimization of ([37]) releasing the 
constraint (|36j). Indeed, computations show that ( pif ) is 
true from t = to t* where the shape of the vortex is 
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determined. Under the assumption that \x and the cur- 
vature do not vary too much along the curve, we derive 
an equation for the minimum 7: 



2a 2 r{t) 2a 2 fL 



p TF {t) (1 + a 2 ) 
2(3 2 z(t) 2a 2 n 



p TF (t) (1 + a 2 ) 



p TF r(t)z(t), 



p TF r(t)r(t). 



Thus, we solve this system with initial conditions r(0) = 
r , r = 0, z(0) = 0, C(0)p TF (0)i(0) = 1. 

We let r$ vary in order to find the minimizing solution. 
We have drawn the vortex line for the minimizing solu- 
tion for some values of £1 in Figure |l} We find that indeed 
the vortex line is bending for a range of f2. The bent vor- 
tex starts to exist near the boundary of the ellipse, that is 
y = v fp Q /a, z = for Q = 21.2, that is &o/to x = 0.368. 
As VL increases, the value of r decreases: r = 0.03 for 
= 21.8, r = 2.9 10~ 4 for Q = 25.8, r = 10~ 6 for 
il = 33.1. As fi increases, r becomes smaller, the bent 
vortex gets very close to the straight vortex. The shape 
of the vortex lines are similar to those obtained in |l!| 
using the full Gross Pitaevskii energy. 




FIG. 1: The vortex line for various values of fl in the z — y 
plane: Q = 21.8 (straight line), Q = 25.8 (dotted line), Q = 
33.1 (dashed line). 

We plot the energy of the straight vortex line and the 
bent vortex vs O in Figure ||. One can observe that for 
Q c = 21.8, that is Q, c /uj x = 0.38 in the initial units, the 
energy of the bent vortex starts to be negative (that is 
below the energy of a solution without vortex), while the 
energy of a straight vortex line is positive. For Q = 33.1, 
the energy of the bent vortex and of a straight vortex 
line become equal. These results are consistent with the 
ones in O] . They obtain the same value of fl c for which 
the bent vortex has a negative energy. 

Let us point out that the bent vortex is a minimizcr 
even if the cross section is a disc. Nevertheless, when e 
is fixed, if (3 gets too big, the straight vortex becomes 
the minimizer, which is the case for [3 = 1. Our analysis 
could give the critical value of (3 above which the vortex 
line should be straight. 



22 24 26 28 



32 34 



FIG. 2: The energy vs. f2 curves for the solution with a 
straight vortex (solid line) and a bent vortex (dotted line) . 



We believe that our analysis justifies why in the condi- 
tions of the ENS experiment, when a vortex is nucleated, 
the contrast is not 100% : indeed, a bent vortex has a 
lower energy than a straight vortex. Nevertheless, the 
velocity of nucleation in the experiment is higher than 
our critical angular velocity Q c : we compute the ther- 
modynamical critical velocity whereas the velocity of nu- 
cleation is likely to be closer to the velocity where the 
vortex free solution loses local stability. Once the first 
vortex is obtained experimentally, if is decreased, the 
bent vortex is likely to exist down to f2 . 

When there are several vortices, we were not able to 
find numerically the shape of the vortex lines which min- 
imize ( |34| ) but we believe that our simplified energy is a 
good description of the experiments ^] and the numerics 
]12[ and we hope that it can be easier to handle than the 
full Gross Pitaevskii energy. 



IV. CONCLUSION 

We have obtained a simplified expression ( p8| ) of the 
energy of a minimizing solution of the Gross Pitaevskii 
energy with a vortex line 7 and (|34|) for n vortex lines ji . 
This expression depends on the shape of the vortex line. 
This has allowed us to draw the vortex line for the min- 
imizing solution and compute its energy. We have seen 
that there is a range of rotationnal velocities for which a 
bent vortex line has a lower energy than a straight vortex 
and a vortex free solution. These computations on the 
simplified expression of the energy are in agreement with 
the computations on the full energy |ll], [T^] . 
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